/*This script generates results at the county level data found in the appendices A, B, and C;
Figure A2
Table B1
Table C1
Table C3
Table C4
Table C5
Table C6
Table C7
Table C8
Table C9

Note: Information to build Table C3 is in output of Table 3 of the main text .do file, County_main.do
*/


cd "$data"
set more off

*Bring in the quarter data
use US_county_2010_Turbines_Irrigation_JAERE.dta, replace

*Globals drawn on
global xcontrols trans_length_a pop_density elev elev_std soil stream_3_15 tmean ppt lat_county lon_county
global xwind ave_wind_class max_wind_class

global sumstat share_og share_og_neg20  ave_wind_class max_wind_class cpis irrigated_a crop_a farm_a t_cap_a first_year ave_year trans_length_a pop_density elev elev_std soil stream_3_15 tmean ppt pback sback

/*Figure A2*/
histogram pback18 if share_og_pos100!=. & pback<10000, scheme(s1color) frequency 
			graph export  "$results/figures/figA2_psetback_dist.pdf", replace


/*Table B1*/

foreach stat in count mean sd min max {
	estpost tabstat $sumstat if samp==1, columns(statistics) statistics(`stat') 
	estadd matrix r=e(`stat')
	eststo tab_count_des_`stat'
	}
	
	
esttab tab_count_des_* ,  label main(r)  nostar  mtitles("Count" "Mean" "St. Dev" "Min" "Max") nonote noobs
esttab tab_count_des_*  using "$results/tables/tabB1_des_county.csv", replace  label main(r)  nostar  mtitles("Count" "Mean" "St. Dev" "Min" "Max") nonote noobs

/*Table C1*/
local i=1

foreach Y of varlist ave_wind_class max_wind_class irrigated_a cpis t_cap_a first_year {

	quietly reghdfe `Y' share_og share_og_neg $xcontrols if share_og_pos100!=., absorb(statefp) cl(cluster_km)
	qui sum `Y' if e(sample)==1
	estadd scalar mean = r(mean)
	estadd scalar Nclus =  e(N_clust)
	
	estimates store share_`i'
	
	local i=`i'+1
	}


esttab share_* using "$results/tables/tabC1_County_Og_share.csv",  label drop(lat_county lon_county _cons) scalar("mean Mean DV" "df_a Fixed Effect Count" "Nclus N Clusters") star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) replace

/*Table C3*/
*hazard attempt.
gen turbine=t_cap_a>0
gen fail=first_year
replace fail=2021 if fail==. 

stset fail, failure(turbine)


local i = 1
foreach irr in "irrigated_a" "cpis" "irrigated_a cpis" {
stcox  `irr' farm_a crop_a $xwind $xcontrols i.statefp if share_og_pos100!=.,  cl(cluster_km)
   		estadd scalar Nclus =  e(N_clust)
   		estadd scalar Nfail =  e(N_fail)
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)
	estimates store county_haz_`i'

	local i = `i' + 1


	
	}

esttab county_haz_* using "$results/tables/tabC3_County_haz.csv",   label  keep( irrigated_a cpis) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  se(a3) b(a3) scalars("r2_p Pseudo R-squared" "Nfail No. Failures" "df_a Fixed Effect Count" "Nclus N Clusters") replace

drop turbine fail

/*Table C4*/
gen constant=1  // the constance variable needs to hard coded in
gen time=1 //time argument necessary to run, but still a cross-sectional regression

tab statefp, gen(spanel) // factor variables are not permitted in the code, this creates state fixed effect variable.

global dist 150 // This sets the maximum distance for the spatial error in kilometers, equivalent to 93.2 miles 

	local i = 1
foreach irr in "irrigated_a" "cpis" "irrigated_a cpis" {
	   qui ols_spatial_HAC t_cap_a `irr' farm_a crop_a $xwind $xcontrols constant spanel2-spanel9 if share_og_pos100!=., lat(intptlat10) lon(intptlon10) timevar(time) panelvar(statefp) distcutoff($dist) // 
		qui sum t_cap_a if e(sample)==1
		estadd scalar mean = r(mean)
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)
		estimates store county_spatial_t_cap_`i'


	qui ols_spatial_HAC first_year `irr' farm_a crop_a $xwind $xcontrols constant spanel2-spanel9 if share_og_pos100!=. & t_cap_a>0, lat(intptlat10) lon(intptlon10) timevar(time) panelvar(statefp) distcutoff($dist) // 
		qui sum first_year if e(sample)==1
		estadd scalar mean = r(mean)
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)
		estimates store county_spatial_fy_`i'
	
	local i = `i' + 1
	
		
	}
	
 	esttab county_spatial_t_cap* county_spatial_fy* using "$results/tables/tabC4_County_Irr_SpatialSE.csv",   label  keep( irrigated_a cpis) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)   se(a3) b(a3) scalars( "mean Mean DV" "df_a Fixed Effect Count") replace	
	

	drop constant time spanel*
	
/*Table C5*/
local s = 0
foreach irr in "irrigated_a" "cpis" "irrigated_a cpis" {
	local s = `s' + 1
	local i = 1
	foreach rest in "share_og_pos100>0" "share_og_pos100!=." "share_og_pos100!=. & statefp!=48" "share_og>0.5 & share_og!=." {
	
		qui tobit t_cap_a `irr' farm_a crop_a $xwind $xcontrols i.statefp if `rest', ll(0) vce(cl cluster_km)
			estadd scalar Nclus =  e(N_clust)
			qui sum t_cap_a if e(sample)==1
			estadd scalar mean = r(mean)
			qui tab statefp if e(sample)==1
			estadd scalar df_a = r(r)
		estimates store sample_t_cap_`s'_`i'
		
		qui reg first_year `irr' farm_a crop_a $xwind $xcontrols i.statefp if `rest' & t_cap_a>0,  vce(cl cluster_km)
			estadd scalar Nclus =  e(N_clust)
			qui sum t_cap_a if e(sample)==1
			estadd scalar mean = r(mean)
			qui tab statefp if e(sample)==1
			estadd scalar df_a = r(r)
		estimates store sample_first_year_`s'_`i'

		local i = `i' + 1
	
		}
	}
	
estimates drop sample_t_cap_2_1 sample_t_cap_3_1  sample_first_year_2_1 sample_first_year_3_1 // this drops the CPIS estimates that are not complete due to incomple CPIS data coverage

*Each panel (A,B, and C) is exported separately*
	
	esttab sample_t_cap_1_* sample_first_year_1_* using "$results/tables/tabC5a_county_irr_sampe.csv",  label  keep( irrigated_a cpis ) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("r2_p Pseudo R-squared" "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count" "Nclus N Clusters") replace

	esttab sample_t_cap_2_* sample_first_year_2_* using "$results/tables/tabC5b_county_irr_sampe.csv",    label  keep( irrigated_a cpis ) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("r2_p Pseudo R-squared" "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count" "Nclus N Clusters") replace
	
	esttab sample_t_cap_3_* sample_first_year_3_* using "$results/tables/tabC5c_county_irr_sampe.csv",    label  keep( irrigated_a cpis ) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("r2_p Pseudo R-squared" "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count" "Nclus N Clusters") replace
	
/*Table C6*/

local s = 0
foreach irr in "irrigated_a" "cpis" "irrigated_a cpis" {
	local s = `s' + 1
	local i = 1
	foreach combo in "i.statefp" "farm_a crop_a i.statefp" "farm_a crop_a $xwind i.statefp" "farm_a crop_a $xwind $xcontrols i.statefp" "$xwind $xcontrols i.statefp " {
	
		qui tobit t_cap_a `irr' `combo' if share_og_pos100!=., ll(0) vce(cl cluster_km)
			estadd scalar Nclus =  e(N_clust)
			qui sum t_cap_a if e(sample)==1
			estadd scalar mean = r(mean)
			qui tab statefp if e(sample)==1
			estadd scalar df_a = r(r)	
		
		estimates store cov_t_cap_`s'_`i'
	
	
		qui reg first_year `irr' `combo' if share_og_pos100!=. & t_cap_a>0,  vce(cl cluster_km)
			estadd scalar Nclus =  e(N_clust)
			qui sum t_cap_a if e(sample)==1
			estadd scalar mean = r(mean)
			qui tab statefp if e(sample)==1
			estadd scalar df_a = r(r)
			
		estimates store cov_first_`s'_`i'
		
		local i = `i' + 1
		
		
	
		}
	}
	
*Each panel (A,B, and C) is exported separately*
	
	esttab cov_t_cap_1* cov_first_1* using "$results/tables/tabC6a_count_irr_cov.csv",  label  keep( irrigated_a cpis ) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("r2_p Pseudo R-squared"  "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count"  "Nclus N Clusters") replace
	
	esttab cov_t_cap_2* cov_first_2*  using "$results/tables/tabC6b_count_irr_cov.csv",   label  keep( irrigated_a cpis ) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("r2_p Pseudo R-squared"  "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count"  "Nclus N Clusters") replace
	
		esttab cov_t_cap_3* cov_first_3*  using "$results/tables/tabC6c_count_irr_cov.csv",   label  keep( irrigated_a cpis ) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalars("r2_p Pseudo R-squared"  "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count"  "Nclus N Clusters") replace 
		
		
/*Table C7*/
*setbacks as outcomes*
qui tobit pback irrigated_a cpis farm_a crop_a $xwind $xcontrols i.statefp if share_og_pos100!=., ll(0) vce(cl cluster_km)
	estadd scalar Nclus =  e(N_clust)
	qui sum pback if e(sample)==1
	estadd scalar mean = r(mean)
	
	qui tab statefp if e(sample)==1
	estadd scalar df_a = r(r)	
		
	estimates store pback

qui tobit pback irrigated_a cpis farm_a crop_a $xwind $xcontrols i.statefp if share_og_pos100!=. & pback<4000, ll(0) vce(cl cluster_km)
	estadd scalar Nclus =  e(N_clust)
	qui sum pback if e(sample)==1
	estadd scalar mean = r(mean)
	
	qui tab statefp if e(sample)==1
	estadd scalar df_a = r(r)	
		
	estimates store pback_out

qui tobit sback irrigated_a cpis farm_a crop_a $xwind $xcontrols i.statefp if share_og_pos100!=. , ll(0) vce(cl cluster_km)
	estadd scalar Nclus =  e(N_clust)
	qui sum pback if e(sample)==1
	estadd scalar mean = r(mean)
	
	qui tab statefp if e(sample)==1
	estadd scalar df_a = r(r)		
	
	estimates store sback
	
*setbacks as covariates*	
local i = 1
foreach combo in "pback sback" "irrigated_a  pback sback" "cpis pback sback" "irrigated_a cpis pback sback"  {
	qui tobit t_cap_a `combo' $xwind $xcontrols i.statefp if share_og_pos100!=., ll(0) vce(cl cluster_km)
		estadd scalar Nclus =  e(N_clust)
		qui sum t_cap_a if e(sample)==1
		estadd scalar mean = r(mean)	
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)	
	estimates store back_cap_`i'
		
		local i = `i' + 1
	}

**remove the controls of setbacks but maintain the sample to observe implied Omitted Variable Bias, not included in preinted table, but referenced in the text
foreach combo in  "irrigated_a" "cpis" "irrigated_a cpis"  {
	qui tobit t_cap_a `combo' $xwind $xcontrols i.statefp if share_og_pos100!=. & pback!=., ll(0) vce(cl cluster_km)
		estadd scalar Nclus =  e(N_clust)
		qui sum t_cap_a if e(sample)==1
		estadd scalar mean = r(mean)	
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)	
	estimates store back_cap_`i'
		
		local i = `i' + 1
	}
		

esttab *back* using "$results/tables/tabC7_county_setback.csv",  label  keep( irrigated_a cpis pback18 sback18) order(irrigated_a cpis pback18 sback18 ) star(* 0.1 ** 0.05 *** 0.01)  se(a3) b(a3) scalars("r2_p Pseudo R-squared"  "mean Mean DV"  "N_lc Censored Observations" "df_a Fixed Effect Count"  "Nclus N Clusters" ) replace


/*Table C8*/
*changes from 1959 to 2007*
local i = 1 
foreach v of varlist perdeltafarmvalue_a perdeltacrop_val_a perdeltaave perdeltattenant perdeltafarm perdeltacrop perdeltafarmintens  {

	
	reg `v' deltacpis deltairr $xwind $xcontrols i.statefp if share_og_pos100!=. ,  cl(cluster_km)
		qui sum `v' if e(sample)==1
		estadd scalar mean = r(mean)
		estadd scalar Nclus =  e(N_clust)
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)	
	estimates store change_`i'
	
	local i = `i' +1
	}

*levels in 2007
local i = 1 
foreach v of varlist lnfarmvalue_a lncropvalue_a lnave_farm ttenant_a  farm_a crop_a farmintens {
	
	
	reg `v' irrigated_a cpis $xwind $xcontrols i.statefp if share_og_pos100!=. ,  cl(cluster_km)
		qui sum `v' if e(sample)==1
		estadd scalar mean = r(mean)
		estadd scalar Nclus =  e(N_clust)
		qui tab statefp if e(sample)==1
		estadd scalar df_a = r(r)	
	estimates store levels_`i'
	
	local i = `i' +1
	}
	
*Each panel (A and B) is exported separately*

	esttab change_* using "$results/tables/tabC8a_County_farmstructure.csv",  label keep(deltacpis deltairrigated) order(deltairrigated deltacpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalar("mean Mean DV" "Nclus N Clusters" "df_a Fixed Effect Count" ) replace
			
	esttab levels_* using "$results/tables/tabC8b_County_farmstructure.csv",  label keep(irrigated_a cpis) order(irrigated_a cpis) star(* 0.1 ** 0.05 *** 0.01)  ar2 se(a3) b(a3) scalar("mean Mean DV" "Nclus N Clusters" "df_a Fixed Effect Count") replace
	
/*Table C9*/
local s = 0
foreach irr in "irrigated_a" "cpis" "irrigated_a cpis" {
	local s = `s' + 1
	local i = 1
	foreach combo in "lnfarmvalue_a lncropvalue_a" "lnave_farm farmintens farm_a crop_a ttenant_a" "lnfarmvalue_a lncropvalue_a lnave_farm farmintens farm_a crop_a ttenant_a "  {
	
		qui tobit t_cap_a `irr' `combo' $xwind $xcontrols  i.statefp if share_og_pos100!=., ll(0) vce(cl cluster_km)
			estadd scalar Nclus =  e(N_clust)
			qui sum t_cap_a if e(sample)==1
			estadd scalar mean = r(mean)
			qui tab statefp if e(sample)==1
			estadd scalar df_a = r(r)	
		
		estimates store t_cap_struct_`s'_`i'
	
		
		local i = `i' + 1
		
			
		}
	}
	
	esttab t_cap_struct_1* t_cap_struct_2* t_cap_struct_3* using "$results/tables/tabC9_County_farmstructure_wind.csv",   label  keep(irrigated_a cpis lnfarmvalue_a lncropvalue_a lnave_farm farmintens farm_a crop_a ttenant_a) order(irrigated_a cpis lnfarmvalue_a lncropvalue_a lnave_farm ttenant_a  farm_a crop_a farmintens) star(* 0.1 ** 0.05 *** 0.01)  se(a3) b(a3) scalars("r2_p Pseudo R-squared"  "mean Mean DV" "N_lc Censored Observations" "df_a Fixed Effect Count"  "Nclus N Clusters") replace
	
********THIS COMPLETES THE COUNTY LEVEL OUTPUT***************
